/ AJ ~ 3 y - C/C 

COMPUTATION OF THERMOCHEMICAL NONEQUILIBRIUM FLOWS 
AROUND A SIMPLE AND A DOUBLE ELLIPSE 

Tahir Gokgen 

Eloret Institute, 3788 Fabian Way, Palo Alto, California 94303 
Work performed at NASA Ames Research Center, Moffett Field, CA 94035 
under Cooperative Agreement NCC2-420 


I. INTRODUCTION 

The nonequilibrium viscous reactive flows over a simple and a double ellipse at 30° an- 
gle of attack have been computed. The geometry and the free stream conditions are given 
by INRIA/GAMNI/SMAI workshop test cases 6.2-2 and 6.2-4. 1 The governing Navier- 
Stokes equations coupled with thermochemical nonequilibrium processes are solved numer- 
ically using a fully coupled, implicit, finite volume technique with a dynamically adaptive 
grid. The nonequilibrium gas model and the numerical method used in the calculations will 
be briefly described in the following sections, while the details can be found in Refs. 2-3. 

II. FORMULATION 

1. Nonequilibrium Gas Model For Air 

The present nonequilibrium gas model for air consists of the following five chronical 
species: three diatomic species, molecular nitrogen IV 2 , molecular oxygen O 2 , nitric oxide 
NO] and two atomic species, nitrogen N and oxygen O. 

Since the diatomic species have internal structure, excitations of the internal degrees of 
freedom, vibrational and rotational, are considered as thermal nonequilibrium processes. 
It is assumed that the local thermodynamic equilibrium within each class of degrees of 
freedom is maintained. 4 ’ 5 These classes are assumed to be the class of translational degrees 
of freedom for all species and the classes of vibrational and rotational degrees of freedom for 
the diatomic species. As a result of the assumed thermal equilibrium within each class, the 
local thermodynamic state of the gas can be described by three temperatures corresponding 
to each energy mode, i. e., translational temperature T, vibrational temperature T v , and 
rotational temperature T r . Equilibrium between these three energy modes is assumed to 
be governed by rate processes. 

The implicit assumption in this three temperature model is that the equilibration time 
scale for each class of the internal mode is small compared to the flow time scale. As a 
result, each internal energy mode is distributed to its equilibrium distribution instanta- 
neously, and can be described by a temperature. 

2. Governing Equations 

The governing equations are basically the continuum Navier-Stokes equations of gas 
dynamics. They are augmented to account for the nonequilibrium processes and are briefly 
presented herein. As a result of the nonequilibrium processes considered, six additional 
partial differential equations are added to the conventional Navier-Stokes equations. The 
equation of conservation of mass is replaced by five spec ies eq uations due to the ch emical 
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nonequilibrium , and two additional energy conservation equations are considered for the 
vibrational and the rotational energy due to the thermal nonequilibrium. 

For simplicity the governing equations are presented in Cartesian tensor notation. The 
species continuity equation, which states the conservation of mass for species i in the gas 
mixture, is given by 

^ + -S-(piUk+jik)=Wi, (1) 

dt oxk 

and the familiar mass-averaged momentum equation, which states the conservation of 
momentum for the gas mixture, is 

?££- + JL(pu j u t + T it +p6 jk )=0. ( 2 ) 

Ot UXk 

Two additional energy equations are introduced to account for the thermal nonequilibrium 
processes, i. e., the relaxation of vibrational and rotational degrees of freedom. Hence, the 
conservation of vibrational energy within diatomic species is expressed by 

^ + -^—{u k E v +?„*) =w v , (3) 

dt oxk 

and similarly, the conservation of rotational energy is 

^ Er + qrk) = ™ r- ( 4 ) 

dt uxk 

Finally, the conservation of total energy is stated by 

+ -x — [uk(E + p) + qk + UjTjk] = 0. (5) 

Ot OXk 

The equation of state or pressure of the gas mixture is given by the Dalton s Law of 
Partial Pressures, 

ns 

p=^ Pj =pfiT. (6) 

t=l 

The total energy per unit volume, E, is made up of the various modes of internal energy, 
i. e., translational, vibrational and rotational modes, the kinetic energy of the flow, and the 
energy arising from the formation of species. These terms, written in the order mentioned, 
give the following expression for the total energy 

1 nS 

E = pC v T + E v + E r + -p{u 2 + v 2 ) + Y,Pi h i ° • ( 7 ) 

1 «=i 

For the diatomic species: the vibrational energy is calculated assuming a harmonic 
oscillator model with a cut off energy level corresponding to the dissociation energy ( trun- 
cated harmonic oscillator ) 6 ’ 7 ; and the rotational energy is calculated assuming a rigid 


dumbell model, 
expressed by 


The vibrational and rotational energies in diatomic species i are then 


_ R , 
Evi = Pi tt( 


&di 


Eri ~ pi Mi Trj 


— 1 e&dilT v _ \ 

R 


), 


( 8 ) 


and the total vibrational and rotational energy in the gas is the sum of the energies in the 
diatomic species. 

Since the source terms associated with the chemical nonequilibrium processes are 
well known (e.g., see Park. 9-11 ), it will not be repeated here. The rate constants for 
chemical reactions are taken from Park. 9 Similarly, in the vibrational and rotational energy 
equations, Eqs. (3) and (4), the source terms w v and w r involve several collisional energy 
exchange mechanisms. 3 ’ 8-11 

The source term in the vibrational energy equation for molecular species i is given by 


Wvi = 


EU ~ E v 


— Wfi Ei(T,T v ) + i obi Gi(T). 


( 9 ) 


The first term on the right hand side of Eq. (9) represents the rate of vibration-translation 
energy exchange for the Landau- Teller model of the harmonic oscillator. Other terms 
are due to the vibration-dissociation coupling (see Ref. 6-7). In the vibration-dissociation 
model, E X (T, T„) is the average energy lost from vibration to a single dissociation and 
Gi(T) is the average energy gained from a recombination. The expressions for E<(T, T v ) 
and Gi(T) are taken from Marrone and Treanor. 7 

The source term in the rotational energy equation is similar to the term in the vi- 
brational energy equation. In the present work, only the rotation-translation exchange 
mechanism, 

Zj i* Li 1 

( 10 ) 


W r 


EU ~ E r 


is considered. 

The diffusive fluxes are assumed to be linearly related to the gradients of the associated 
flow quantities. The Navier-Stokes diffusion of momentum relates the shear stress tensor 
to the velocity gradient tensor, i. e., the Newtonian fluid assumption; Fourier’s law of heat 
conduction relates the heat flux vector to the associated temperature gradient; and Fick’s 
law of diffusion relates the diffusive mass fluxes to gradients of species concentration. 

For the viscosity of species i the following expression is used: 


ln/ii = (>l t TnT + 2?i)lnT + Ci . 


( 11 ) 


The coefficients A;, S,, and Ci for species i are taken from Moss 12 . The thermal conduc- 
tivity of translational energy for species i is calculated using 


Ki — — Hi C V(i . 


( 12 ) 



Once the viscosity and thermal conductivity of each chemical species are specified, the 
viscosity and thermal conductivity of the gas mixture are computed using the semiempirical 
Wilke’s mixture rule. 13 The thermal conductivities of vibrational and rotational energies, 
analogous to Eucken relation, are calculated by 


nm 

i=l 

nm 

Kr 1=1 /*i* 

i—l 


(13) 


The diffusion coefficient V is specified assuming constant Schmidt number for all 
species, 

&=-^r = 0.5. (14) 

pV 

The vibrational and rotational relaxation times are specified such that the collision 
numbers of the vibrational and rotational relaxation, Z v and Z r , are 50 and 5, respectively. 
The mean collision time is estimated using the mean free path A and the mean speed c 


X 

Tcoll. — - > 

C 

and the vibrational and rotational relaxation times are then calculated by 


(15) 


T„ = Z v T co U . = 50 T C 0 U, 
T r = Z r T co \i . = 5 T co jj. 


(16) 


The above specifications are made such that the present nonequilibrium air model is com- 
patible with the five species DSMC air model of Moss et a/. 14 


m. METHODOLOGY 

1. Numerical Method 

The numerical method is based on that presented by MacCormack 15 for solving the 
two dimensional perfect gas Navier-Stokes equations. This method has the following useful 
features: it is implicit, which allows one to take large time steps to reach the steady 
state solution; it uses a flux vector splitting similar to that of Steger and Warming 16 to 
mimic convective transport and/or to maintain numerical stability; and it uses the Gauss- 
Seidel Line Relaxation technique with the block tridiagonal matrix inversion. Candler 
and MacCormack 17-19 applied this method to solve the two dimensional Navier-Stokes 
equations with thermochemical nonequilibrium processes. In the present work, the method 
is extended to include a shock capturing technique with dynamically adaptive grids. 

2. Adaptive Grid Algorithm 

The adaptive grid strategy used here follows the work presented by Eiseman 20 and 
Abolhassani et al. 21 Grid motion is established directly in one physical dimension (normal 
to the wall) using a weight function in such a way that the weight function and some 


measure of the mesh size would be held approximately constant for the grid. Large values 
of the weight cause small mesh spacing and produce the desired clustering of the points 
while small values of the weight produce the opposite effect, i.e., large mesh spacing. The 
useful form of the weight function is given by the linear combination 20 

m 

w = l + Y J C k M k . (17) 

k= 1 

As to the choice of the weight functions M*, the gradients of the flow variables seem 
to be a natural choice. This choice produces good results as long as the chosen flow 
variable is monotonic in the direction that the grid is being adapted. Otherwise, a simple 
gradient produces coarse mesh spacing where the chosen flow variable has an extremum. 
This disadvantage can be avoided by incorporating the second derivative into the weight 
function. In the present work, five functions for axe used. These are the gradients of 
the pressure, tangential velocity, species mass fraction, temperature, and kinetic energy of 
the flow. These functions proved to be satisfactory for a wide range of flow conditions. 

IV. RESULTS 

The results presented here axe for the viscous reactive flows over a simple and a 
double ellipse at 30° angle of attack. The simple/double ellipse geometry is provided by 
INRIA/GAMNI/SMAI workshop test cases 6.2. 1 For all calculations, a freestream Mach 
number Moo = 25, unit Reynolds number Rtoo/m = 2.2 x 10 4 , and temperature and 
pressure of 75 km standard atmosphere (T = 205.3°AT and p = 2.52 Pa) are prescribed. 
The wall is assumed to be at 1500°# and noncatalytic to chemical reactions. 

The adaptive grid used for the simple ellipse case (test case 6.2-2) is depicted in Fig. 1. 
As seen from Fig. 1, the grid around the ellipse is adapted to the flow gradients of the 
bow shock wave and the viscous boundary layer near the wall. In Figs. 2-8, the flowfield 
over the ellipse is represented in terms of contour plots. Fig. 2 shows the Mach number 
contours. It should be mentioned that the Mach number is based on the frozen speed of 
sound. Since the rotational degrees of freedom are treated to be nonequilibrium with the 
translational degrees in the gas, the speed of sound is larger than that of one (or two) 
temperature models in which translational and rotational degrees are assumed to be in 
equilibrium. Therefore, the computed Mach number is lower than that of one (or two) 
temperature models. In Figs. 3-5, translational, rotational and vibrational temperature 
fields are shown, respectively. It appears, from Figs. 3 and 4, that the translational and 
rotational temperature fields are identical. For the flow conditions prescribed in this test, 
the rotational relaxation occurs very fast behind the shock wave. In fact, everywhere in 
the flowfield, except for a few points right behind the shock, the rotational temperature is 
in equilibrium with the translational temperature. This justifies the assumption that one 
temperature can characterize translational and rotational modes of the internal energy 
at low altitudes. However, for more rarefied cases this assumption breaks down, and 
the rotational relaxation needs to be considered. As far as the vibrational relaxation is 
concerned, a close examination of Fig. 4 and 5 shows that the vibrational temperature 
in most of the flowfield is different from the translational temperature. The previous 
temperature figures show the extent of thermal nonequilibrium in the flowfield. Figs. 6 



and 7 show the extent of chemistry in the flowfield in terms of atomic nitrogen and atomic 
oxygen mass fraction contours. Fig. 6 shows that considerable amount of N is produced 
due to dissociation of N 2 near the nose region and on the compression (windward) side of 
ellipse and to a smaller extent on the expansion (leeward) side. It is observed from Fig. 7 
that almost all 0 2 is dissociated right behind the shock wave and recombined in the cool 
boundary layer on the compression side. 

Most of the comments made on the simple ellipse calculations equally apply to the 
double ellipse calculations. Fig. 8 shows the grid used for the double ellipse calculations. 
Again, the converged grid is adapted to the flow gradients of the bow shock wave and the 
viscous bo un dary layer near the wall. The computed flowfield around the double ellipse is 
represented by the contour plots Figs. 9-13, and the surface quantities Figs. 14-18. Fig. 9 
shows the pressure coefficient contours. Although the detached bow shock wave is obvious 
from the contours, the canopy shock is not observed because of the increment used in 
the contours of C p . Figs. 10-11 represents the translational and vibrational temperature 
contours and the extent of thermal nonequilibrium in the flow, and Figs. 12-13 give the 
contours of atomic nitrogen and oxygen. The surface quantities, the pressure coefficient 
C p , the skin friction coefficient Cf and the Stanton number St, are shown in Figs. 14-16, 
respectively. Finally, distributions of atomic species N and O at the wall are shown in 
Figs. 17 and 18. 

V. CONCLUDING REMARKS 

The thermochemical nonequilibrium viscous flows are solved numerically using a fully 
coupled, implicit, finite volume technique with a dynamically adaptive grid. The adaptive 
grid is necessary because there are regions of the flowfield where a predetermined grid fails 
to resolve the flow characteristics such as near shock waves and the viscous boundary layer. 

The present air model includes thermal as well as chemical nonequilibrium effects. 

It is shown that most of the flowfield is in vibrational nonequilibrium, which may show 
significant differences with respect to flow chemistry compared to one temperature model 
suggested in the workshop problem. However, it is clear that for the computed test cases 
the rotational relaxation is not important. 
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Fig. 1. Adaptive grid for reactive flow over a simple ellipse (60 x 45). 
Test Case 6. 2-2: a =30°, = 25, Re^/m = 2.2 x 10 4 
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Fig. 2. Mach number contours over a simple ellipse. 

Test Case 6. 2-2: a = 30°, = 25, Re^/m = 2.2 x 10 4 












Fig. 4. Rotational temperature contours over a simple ellipse. 

Test Case 6. 2-2: a - 30°, = 25, Re^/m = 2.2 x 10 4 




Fig. 5. Vibrational temperature contours over a simple ellipse. 

Test Case 6.2-2: a = 30°, Moo = 25, Re^/m = 2.2 x 10 4 
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Fig. 6. Atomic nitrogen N contours over a simple ellipse. 

Test Case 6. 2-2: a = 30°, = 25, Re^/m = 2.2 x 10 4 








Fig. 7. Atomic oxygen O contours over a simple ellipse. 

Test Case 6. 2-2: a = 30° , = 25, Re^/m = 2.2 x 10 4 




Fig. 8. Adaptive grid for reactive flow over a double ellipse (90 x 60). 
Test Case 6. 2-4: a = 30°, M 0 0 = 25, Re oo/m = 2.2 x 10 4 




Fig. 9. Pressure coefficient contours over a double ellipse. 

Test Case 6. 2-4: a = 30°, = 25, Re^/m = 2.2 x 10 4 






Fig. 10. Translational temperature contours over a double ellipse. 

Test Case 6. 2-4: a = 30°, = 25, Re^/m = 2.2 x 10 4 





Fig. 11. Vibrational temperature contours over a double ellipse. 

Test Case 6. 2-4: a = 30°, = 25, Re^/m = 2.2 x 10 4 





Fig. 12. Atomic nitrogen N contours over a double ellipse. 

Test Case 6. 2-4: a = 30°, = 25, Re^/m = 2.2 x 10 4 







v_y 


0.7 



Fig. 13. Atomic oxygen O contours over a double ellipse. 

Test Case 6. 2-4: a = 30°, = 25, Re^/m = 2.2 x 10 4 







Fig. 14. Distribution of pressure coefficient at the wall of the double ellipse. 
Test Case 6. 2-4: a = 30°, = 25, Re^/m = 2.2 x 10 4 
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Fig. 15. Distribution of skin friction coefficient for the double ellipse. 
Test Case 6. 2-4: a — 30°, M 0 0 = 25, Re oo/m = 2.2 x 10 4 
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Fig. 17. Distribution of N mass fraction at the wall of the double ellipse. 
Test Case 6. 2-4: a = 30°, = 25, Re^/m = 2.2 x 10 4 





